program cubeflu

use basic_var
use mpi_basic_var
use timing
use mpi
#ifdef _ACCEL
use accel_lib
#endif
#ifdef _OPENMP
use omp_lib
#endif

real(myk) :: start_time,end_time,elapsed_time

call start_mpi()

GPU = .false.

#ifdef _ACCEL
print*,'NEW n_rank,NCOREXNODE: ',n_rank,NCOREXNODE
if(mod(n_rank,NCOREXNODE) .lt. NGPUXNODE) then
  call acc_set_device(acc_device_nvidia)
  device_kind = acc_get_device()
  print*,'Selected device kind:  ',device_kind
  num_devices = acc_get_num_devices(acc_device_nvidia)
  print*,'Number of devices available: ',num_devices
  call acc_set_device_num(mod(n_rank,NCOREXNODE),acc_device_nvidia)
  print*,'n_rank: ',n_rank,' tries to set GPU: ',mod(n_rank,NCOREXNODE)
  my_device = acc_get_device_num(acc_device_nvidia)
  print*,'n_rank: ',n_rank,' is using device: ',my_device
  print*,'Set GPU to true for rank: ',n_rank
  GPU = .true.

  print*,'FOR POST-PROCESSING SET GPU TO FALSE IN ANY CASE'
  GPU = .false.
endif
#endif
print*,'FOR POST-PROCESSING SET OPENMP TO FALSE IN ANY CASE BUT MODIFY NCOREXNODE TO NGPUXNODE+1 (i.e.=1)'
if(openmp_activate == 1) then
  NCOREXNODE = NGPUXNODE + 1
endif
openmp_activate = 0
#ifdef _OPENMP
if(.not. GPU) then
   if(openmp_activate == 1) then
     print*,'Setting OpenMP threads :',NCOREXNODE - NGPUXNODE
     call omp_set_num_threads(NCOREXNODE - NGPUXNODE) ; NCOREXNODE = NGPUXNODE + 1
   else ; call omp_set_num_threads(1) ; endif
endif
#else
if(openmp_activate == 1) then
  print*,'Error: please compile using OpenMP or disable openmp_actiavte in mpi_gpu.dat' ; call end_mpi() ; STOP
endif
#endif

call post_init()
call size_init()
call decompose_mpi()
call allocate_storage()
print*,'MPI and Memory allocated succesfully'

if(MPI_STAT > 0)  call input_grid()
call input_grid_tot()

if(GPU) then
!  call update_gpu()
  print*,'GPU activated for rank : ',n_rank
endif
call mpi_barrier(mpi_comm_world,ierr)
call time(start_time)

if(MPI_STAT > 0) then
   call init_pos()
   ! First order statistics
   if(n_rank == 0) print*, '__________ I ORDER _____________'
   do post_it=post_itmin,post_itmax,post_itout
      call read_sol_binary_pio()
      call print_screen()
      if(n_rank.eq.0) then  
          print*, '__________________________________________________'
          print*, 'post_it = ',post_it,' t = ',t
      endif 
      if(MPI_STAT == 4) call print_vtk()
      call post_mpi_first()
   
   enddo
   call reduce_mpi_space_and_time_first()
   
   ! Second order statistics
   if(n_rank == 0) print*, '__________ II ORDER _____________'
   do post_it=post_itmin,post_itmax,post_itout
      call init_pos_space_average()

      call read_sol_binary_pio()
      call print_screen()
      if(n_rank.eq.0) then 
          print*, '__________________________________________________'
          print*, 'post_it = ',post_it,' t = ',t
      endif  
      call post_mpi_second()

      call reduce_mpi_space_second()

      if(n_rank .eq. 0) call print_average_tecplot_space(1)
   enddo
   post_it = post_itmax ! to reduce using post_it instead of post_itmax
   call reduce_mpi_space_and_time_second()
   
   ! Print averaged fields
   if(n_rank .eq. 0) call print_average_tecplot_space_and_time(1)
endif

! Full field statistics
if(FULL_STAT > 0) then
   if(n_rank.eq.0) then  
      call init_tot()
      call init_pos()
      if(FFT_STAT > 0) call init_fft()
      if(n_rank == 0) print*, '__________ FULL: I ORDER _____________'
      do post_it=post_itmin,post_itmax,post_itout
!!         call read_sol_binary()
!!         call reconstr_MPI()
         call read_sol_binary_pio_stream()
         print*,'Full field reconstructed'
         print*, '__________________________________________________'
         print*, 'post_it = ',post_it,' t = ',t
         if(FULL_STAT == 2)  call total_print_tecplot()
         if(FULL_STAT == 3)  call total_print_plot3d()
         if(FULL_STAT == 4)  call total_print_vtk()
         if(FULL_STAT == 10) call write_sol_binary_pio_stream()
         call post_full_first()
         if(FFT_STAT > 0) then
           call transform_fft()
           call post_fft()
         endif
      enddo
      call reduce_full_space_and_time_first()
      if(FFT_STAT > 0) then
        call reduce_fft()
        call print_fft()
      endif
      
      do post_it=post_itmin,post_itmax,post_itout
         call init_pos_space_average()
!!         call read_sol_binary()
!!         call reconstr_MPI()
         call read_sol_binary_pio_stream()
         print*,'Full field reconstructed'
         print*, '__________________________________________________'
         print*, 'post_it = ',post_it,' t = ',t
         call post_full_second()

         call reduce_full_space_second()

         call print_average_tecplot_space(2)
      enddo
      post_it = post_itmax ! to reduce using post_it instead of post_itmax
      call reduce_full_space_and_time_second()

      ! Print averaged fields
      call print_average_tecplot_space_and_time(2)
   endif
endif

call time(end_time)
elapsed_time = end_time - start_time
print*,'Elapsed time : ',elapsed_time,' GPU? ',GPU

call end_mpi()

end program cubeflu
